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Abstract 

We extend a non-Tikhonov asymptotic embedding, proposed earlier, for calculation of 
conduction velocity restitution curves in ionic models of cardiac excitability. Conduction 
velocity restitution is the simplest nontrivial spatially extended problem in excitable media, 
and in the case of cardiac tissue it is an important tool for prediction of cardiac arrhythmias 
and fibrillation. An idealized conduction velocity restitution curve requires solving a nonlinear 
eigenvalue problem with periodic boundary conditions, which in the cardiac case is very stiff 
and calls for the use of asymptotic methods. We compare asymptotics of restitution curves in 
four examples, two generic excitable media models, and two ionic cardiac models. The generic 
models include the classical FitzHugh-Nagumo model and its variation by Barkley. They are 
treated with standard singular perturbation techniques. The ionic models include a simplified 
"caricature" of the Noble (1962) model and the Beeler and Renter (1977) model, which lead 
to non-Tikhonov problems where known asymptotic results do not apply. The Caricature 
Noble model is considered with particular care to demonstrate the well-posedness of the 
corresponding boundary-value problem. The developed method for calculation of conduction 
velocity restitution is then applied to the Beeler- Renter model. We discuss new mathematical 
features appearing in cardiac ionic models and possible applications of the developed method. 

Keywords action potential; traveling wave. 
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1 Introduction 



Cardiac excitability models Hodgkin and Huxley's model of the electric properties of the 
giant squid axon [29] was the first to describe in mathematical terms the exclusively biological 
phenomenon of excitability. It started a revolution in science well-worth the Nobel prize it was 
awarded. This achievement has been followed by the development of a long sequence of mathe- 
matical models of heart excitability starting from Noble's works [48, 49]. Due to its importance 
for biomedical applications, particularly for understanding and treatment of cardiac arrhythmias 
caused by pathologies of electrical excitation and propagation, the mathematical modelling di- 
rection has been under intensive development during the last decades, and currently has reached 
clinical applications and industrial scale. It lies in the heart of the ambitious Physiomc project [30] 
which aims at a mathematical description of the physiology of whole organisms. Due to complex- 
ity of the models involved they are mainly used in numerical computations and contribute a 
substantial load on the UK national supercomputer facilities [51]. 

Stiffness of cardiac excitability models The computational complexity of cardiac models 
lies not only in the complexity of the heart as a system, which compared to the brain is rela- 
tively modest, but also in the essential stiffness of cardiac equations. These equations have to 
describe very sharp and fast excitation fronts where some processes happen on the scale of tens 
of microseconds and micrometers, through to tissue and organ level, on the scale of seconds and 
centimeters, thus covering several orders of magnitude. Thus a challenge for applied mathematics 
is how to turn this stiffness from an adversary into an ally. A standard approach is to treat small 
parameters, responsible for such stiffness, using asymptotic rather than numerical methods. For 
the Hodgkin-Huxley model, a simple caricature easily treatable mathematically has been intro- 
duced by FitzHugh [22] and Nagumo et al. [47] , which was based on a modification of the classical 
van der Pol system of equations [59]. Asymptotic analysis of FitzHugh-Nagumo type systems, a 
nice summary of which can be found e.g. in [58], has achieved remarkable success in describing, 
in qualitative terms, many of the phenomena observed in more realistic, experiment based ionic 
models. 

The traditional asymptotic approach The essence of the approach is separation of the 
dynamic variables into "fast and slow" , similar to the classical Tikhonov-Pontryagin scheme [44, 
52, 57], only in a spatially extended context. A typical solution consists of moving, fast and steep 
"fronts" and "backs" of excitation pulses, located near codimension-one manifolds, i.e. points in 
one spatial dimension (ID), lines in two spatial dimensions (2D), and surfaces in three spatial 
dimensions (3D), which are interspersed by smooth and slow intervals. During the fast fronts and 
backs, the slow variables remain almost unchanged. During the slow intervals, the fast variables 
remain very close to their quasi-stationary values determined by the current values of the slow 
variables. The slow pieces are typically of two kinds: with lower and with higher values of the 
transmembrane voltage or a variable that corresponds to it. The lower- voltage, "diastolic" pieces 
are close to or include the "resting state" , representing excitable tissue which was not excited 
for a long time, and the higher-voltage, "systolic" pieces represent the "action potential" phase 
of the excitation. An extra feature of 2D and 3D is the possibility of "wave breaks" , which are 
of particular relevance for cardiac arrhythmias. Such wave breaks are moving codimension-two 
manifolds, i.e. points in 2D and lines in 3D, where the fronts and backs meet. It is essential that 
mathematically, fronts and backs are objects of the same nature, differing only in the direction 
of motion: at the fronts, the systolic phase advances; at the backs, the systolic phase recedes, so 
the wave break is where the interface between the phases is momentarily stalled. This description 
allows even some analytical treatment of the motion of wave breaks in 2D, including steadily 
rotating and meandering spiral waves of excitation [25] . Conceivably, this asymptotic description 
could be also used numerically within an appropriate moving interface methodology. 
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The need for a non-Tikhonov embedding However, sinee models of FitzHugh-Nagumo type 
have been typically postulated rather than derived from realistic ionic models of cardiac excitation, 
the question about their quantitative validity was not usually posed. Successful attempts to 
apply the same singular perturbation technique as developed for systems of FitzHugh-Nagumo 
type, directly to detailed ionic models, have been made, e.g. [38], but this did not turn into a 
mainstream practical approach. We believe that the reason is that systems of FitzHugh-Nagumo 
type are actually quite different, in the asymptotic sense, from detailed ionic cardiac models, as 
they fail completely to describe, even at a qualitative level, some important properties of cardiac 
excitation, such as 

• slow reporarization, 

• slow subthreshold response, 

• fast accommodation, 

• variable peak voltage and 

• front dissipation, 

all of which are experimentally well-established and also successfully reproduced by detailed cardiac 
ionic models [7, 9]. The slow repolarization means that, although cardiac excitation pulses do 
indeed possess steep fronts, they have no steep backs, at least not steep enough compared to the 
steepness of the fronts, anyway. Hence interpretation of wave breaks in 2D and 3D as loci where 
fronts meet backs is inapplicable to cardiac models for absence of backs. Since propagation blocks 
and wave breaks are very important in most applications of mathematical cardiology, there is no 
much hope that the FitzHugh-Nagumo ideology could lead to a practical numerical tool that could 
tame the stiffness of the cardiac equations. 

In a recent series of works [7-9, 11], we have developed an analytical approach to cardiac equa- 
tions based on their special structure, different from the FitzHugh-Nagumo paradigm, and taking 
into account small parameters actually present (sometimes hidden) in the equations, rather than 
trying to force them into the Procrustean bed of a classical scheme. Using the existence of large (or 
small) values of some variables for model reduction and perturbation analysis is a basic technique in 
applied mathematics. One well-known example is the Quasi(Pseudo)-Steady-State approximation 
[26, 54], which reduces the equations of an enzyme reaction to a singular perturbation Tikhonov 
problem. Another prominent example is the wide application of scale separation and model re- 
duction techniques to problems of chemical kinetics and reactive flows. Since in such problems 
the number of reacting species is huge, this is typically done by computational algorithms such as 
Computational Singular Perturbations (CSP), Intrinsic Low-Dimensional Manifolds (ILDM), the 
Grad Moment method and others [24, 33] . It has been shown that these computational techniques 
generate the asymptotic expansion of a slow invariant manifold of a Tikhonov problem [33, 63]. 

However, the application of asymptotic embedding techniques is not restricted to Tikhonov 
problems, nor must it, a priori, lead to such. Indeed, our recent works [7-9, 11] have clearly 
demonstrated that, to achieve a physiologically correct asymptotics in realistic models of cardiac 
excitation, a parameter embedding is needed which involves a large factor in front of individual 
terms, but not the whole, of the right-hand side of some equations (e.g. the /Na term in the trans- 
membrane voltage evolution equation), non-analytical, perhaps even discontinuous, asymptotic 
limit of some right-hand sides (e.g. the I^a gating variables), even though the original system is 
analytical, non-isolated equilibria in the fast subsystem and dynamic variables which change their 
character from fast to slow within one solution (e.g. the transmembrane voltage). 

In particular, we have demonstrated that separate consideration of the fast subsystem describ- 
ing excitation front produces a simple useful criterion of dynamic propagation block in a modern 
cardiac model [55], and the fast and slow subsystems can be successfully matched to describe the 
action potential as a singular limit in a single-cell (OD) variant of a simple cardiac model [9] . The 
next step is to combine the fast and slow description in a spatially extended context. 
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CV restitution curves: a spatiotemporal problem involving fast and slow scales The 

aim of the current work is to make this next step. For this purpose, we have chosen the simplest 
nontrivial spatially extended problem that depends both on the fast and the slow processes: the 
conduction velocity restitution curve. This choice is also motivated by the practical importance of 
restitution curves, which are of two kinds. The action potential duration (APD) restitution curve 
is the dependence of the APD on the duration of the preceding diastolic interval (DI). Nolasco 
and Dahlen [50] noted that in a single-cell setting and with a fixed period of excitation, a slope 
of the APD(DI) curve greater than one indicates instability of the even APD sequence. For this 
reason the restitution curves are considered an important instrument in understanding instabilities 
of excitation waves leading to onset of cardiac arrhythmias. Later studies have demonstrated 
that in a spatially extended context, another important tool is the conduction velocity (CV) 
restitution curve, which describes the dependence of CV on the preceding diastolic interval. The 
CV(DI) dependence together with the APD(DI) dependence and the fact that the overall period 
known in electrophysiology as Basic Cycle Length is given by BCL=APD-|-DI, makes it possible 
to define the CV(BCL) dependence, i.e. relationship between the period of excitation waves and 
their propagation speed, which is also known as the dispersion relation in general wave theory. 
The CV(DI) curve depends on the definition of the boundary between action potential phase 
and the diastolic phase, which for cardiac excitation pulses is arbitrary for lack of sharp backs. 
The CV(BCL) dependence is, on the contrary, free from such arbitrariness and is well defined 
mathematically. So in our study, we shall use this dependence as the restitution curve. 

Types of restitution curves In view of the clinical importance of fibrillation, numerous exper- 
imental, e.g. [12, 14, 23, 50, 61] and numerical e.g. [17, 18, 31, 35, 36, 60] studies are concerned with 
tests of this hypothesis, and with measurements and computation of restitution curves in various 
types of cardiac cells. Measurement and computation of restitution curves are not straightforward. 
A number of different experimental/numerical protocols are in use (see e.g. [53] and references 
therein), which produce different curves and it is not always clear which is the most relevant one 
in a particular case. For instance, in the so called "dynamic" protocol the tissue is paced at a 
given basic cycle length until a periodic regime is established, and the APD, DI and CV of the 
established pulses are recorded. Then the process is repeated with other cycle lengths. Another 
protocol is the "S1-S2" restitution protocol, in which the tissue is paced at a fixed cycle length SI 
until a periodic regime is reached, and is then perturbed by an out-of-sequence stimulus (S2) and 
the response is recorded. The preparation is then paced at a the same SI until steady-state has 
been reached again, and is then perturbed by a different S2. The curve so measured depends on 
the choice of the SI cycle length, and therefore it is not even unique. Although used in electrophys- 
iological practice, these protocols have a number of drawbacks: they contain some arbitrariness 
and thus lead to results which are not unique, they are prone to systematic errors since it not 
easy to distinguish the ultimate periodic regime from transient, they are time consuming and, in 
the case of numerical simulations, computationally expensive since a repeated solution of large 
systems of stiff nonlinear partial differential equations is required. 

The aim of this study In the present study we have chosen to use an idealized definition of 
the "dynamic" restitution protocol, i.e. we consider strictly periodic wave solutions, and study the 
dependence between the BCL and CV of such solutions. This dependence is well defined math- 
ematically via solvability of the corresponding boundary-value problem with periodic boundary 
conditions. This idea is not new but so far it has had only limited application for the following 
two reasons. First, the resulting boundary value problem is typically very stiff, with very steep 
upstroke but slow prolonged plateau and recovery stages of a typical cardiac action potential, and 
so its direct solution requires considerable effort. Secondly, the Tikhonov asymptotic embeddings 
which are typically used to alleviate such scale disparities fail to produce results which are even 
qualitatively correct, as noted above. 

We wish to emphasize that the periodic boundary value problem approach we advocate here 
is applicable both to cardiac models with Tikhonov as well as with non-Tikhonov asymptotic 
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structure and in this work wc illustrate both of these cases. However, in the absence of a rigorous 
theory of the non-Tikhonov case, we make a special effort to investigate whether the resulting 
asymptotic boundary value problem is well-posed. This is not obvious a priori. 



Structure of the paper In section 2 we formulate the periodic boundary-value problem which 
gives a general method of computing CV restitution curves regardless of the asymptotic structure 
of the particular cardiac model. In sections 3, 4 and 5 we apply the method to well-known models 
with Tikhonov asymptotic structure in order to provide simple illustrations. Section 6 is central 
to the article. Here we use a suitably reformulated version of the Noble model of cardiac Purkinje 
fibers [49] to illustrate the non-Tikhonov asymptotic reduction in detail and to investigate whether 
the resulting asymptotic boundary value problem is, indeed, well-posed. In section 7, we calculate 
the full and the asymptotic CV restitution curves of the Beeler-Reuter ventricular model [6] and 
demonstrate a good quantitative agreement. Section 8 provides concluding remarks and suggests 
possible extensions of the work. 

2 Restitution curves: the boundary-value problem formu- 
lation 

A typical voltage-gated model of cardiac excitation and propagation in a one-dimensional, homo- 
geneous and isotropic medium has the form of a reaction-diffusion system. 



where x is the spatial coordinate, t is the time, E is transmembrane voltage of the cardiocytes, the 
functions /; (•) represent individual transmembrane ionic currents, each conducted by a specific type 
of transmembrane channel, the vector y includes a number of "gating" variables controlling the 
permittivity of the ionic channels and the intra- and extracellular concentrations of ions involved, 
and _D > is a "voltage diffusion constant" , depending on the electric capacitance of cardiocytes 
and Ohmic contacts between them. Note that D can be made equal to any positive value by 
rescaling the spatial variable x; we shall choose this scaling so that _D = 1, or related to the small 
parameter when considering asymptotics. This means that the dimensionality of x is that of t^/^. 
Correspondingly, to compare our subsequent results with experimental data, lengths and speeds 
should be scaled up by the factor of D^^^, where D is the value of the voltage diffusion coefficient, 
dependinging on the properties of the given tissue and the direction of wave propagation. 

The number of gating variables, concentrations and the form of the functions //(•) and Fy(-) are 
fitted to reproduce the very latest experimental observations. As experimental methods improve, 
the models evolve to be ever more complicated but the general form of the reaction-diffusion system 
(1) has hardly changed since 1962 when the first cardiac model was published by Noble [49]. A 
relatively recent but by no means ultimate list of cardiac models can be found in the review [15] 
and confirms this assertion. 

CV restitution curves are typically computed by direct numerical simulation of the partial 
differential equations (1) following a particular protocol. As argued above this is computationally 
expensive and time consuming, and prone to systematic errors. A more sound mathematical 
approach is to look for solutions in the form of waves travelling with a constant velocity c > 
and a fixed shape. This is guaranteed by the travelling wave ansatz F{z) = F{x — ct) for the 
dynamical variables F = E,y, where z = x — ct. Equations (1) are then reduced to a system of 
autonomous ordinary differential equations and the CV restitution curve can be found from the 
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periodic boundary value problem 



d^ 



AE 
'd7 



+ c— + ^/i(^,y) =0, 

I 



c^ + Fy(£;,y)=0, 
az 



E{0) = E{cP), ^ 



AE 
Az 



y(0)=y(cP), E{Q)^Ea, 



(2a) 
(2b) 

(2c) 



where P is the temporal period of the waves. The last boundary condition i?(0) = Eq is related 
to the translational invariance of the problem. This condition allows the selection of a single 
solution out of a one-parametric family of solutions differing from each other only in their position 
along the z axis; thus the exact choice of Eq is not essential, as long as it is selected within the 
range of values of E{z). Problem (2) is of order (dim(y) + 2), where dim(y) is the dimension 
of the vector y, and its general solution includes (dim(y) + 2) arbitrary constants. In addition, 
the problem involves two unknown parameters, c and P. On the other hand, it has (dim(y) + 2) 
periodic boundary conditions plus the last "phase" or "pinning" condition required to eliminate 
the translational invariance of the system. Thus, we have (dim(y)-|-4) parameters and (dim(y)-|-3) 
constraints on them, so the solution of the problem should yield, in principle, a one-parameter 
family of solutions. A projection of this family onto the (P, c) plane is the sought after "ideal" 
dynamic CV restitution curve describing the dependence of the wave speed on the wave period. 



3 Outline of the singular perturbation theory of Tikhonov 
excitable systems 

The method outlined in section 2 is applicable to any cardiac model but due to the inherent stiffness 
of cardiac equations solution it is difficult for a numerical study if asymptotics are not exploited. 
Our first illustrations of the method will involve the Barkley model [5] and the FitzHugh-Nagumo 
system [22, 47]. We take advantage of the fact that these models have well-known asymptotic 
structures of a Tikhonov type which has been studied in a number of works, e.g. [19, 37, 39, 41, 
45, 58] and thus they provide simple illustrations and set the context for our main results. 

3.1 Asymptotic reduction of the CV restitution boundary value prob- 
lem for Tikhonov systems 

A typical model with a Tikhonov asymptotic structure has the form 



dE ,^ , d^E 

dy 



FEiE,y)+e—, (3a) 



Fy{E,yl (3b) 

where E is interpreted as the voltage while y is taken to represent all other variables of an ionic 
cardiac model. We assume, for simplicity, that the dynamical variables E and y are scalar fields, 
which is true for the Barkley and the FitzHugh-Nagumo models. The small parameter e <C 1 
specifies the asymptotic structure of the system explicitly by indicating the relative magnitude 
of the various terms in the model. Here and in subsequent asymptotic formulations, the spatial 
scaling is chosen so that the diffusion coefficient is equal to e; the convenience of this choice will 
be evident shortly. In order for equations (3) to have excitable or oscillatory dynamics, certain 
properties of the functions Fe{) and Fy{-) need to be assumed. We shall assume that in a certain 
interval of y values, y EXy ^ (ymin, J/max), the following is true. 
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Al. Equation Fe{E, y) — 0, understood as an equation for E aX a, fixed y, has three roots, namely 
E-{y) < E^,{y) < E^{y). This equation defines the redueed slow manifold of the system (3) 
in the sense of the geometric singular perturbation theory [32, 34]. 

A2. The roots have alternating stability in linear approximation, that is, OeFe iE_(jj),y) < 0, 
OeFe {E^{y),y) > and OeFe {E+{y),y) < 0, where Oe denotes a partial derivative. 

Further, we assume that, in a possibly smaller interval y ^I'y — {y'^i-^, 2/max) ^ ^i/j the following 
is true. 

A3. The slow dynamics for y is growth for the lower root E = E^{y) and decrease for the upper 
root E = E+{y), that is Fy{E_{y),y) > and Fy{E+{y),y) < 0. 

A4. The periodic wave solutions of interest only involve the interval y ^Ty. 

These assumptions are true for the Barkley and the FitzHugh-Nagumo models, and are illus- 
trated in figures 1(a) and 2(a) below. The last assumption A4, unlike the first three, is difficult to 
formulate in a priori terms, and we shall discuss its implications as we obtain the relevant results 
below. 

Due to assumption Al, the sets {E-{y),y)) and {E+{y),y)) are disjoint in the {E,y)-p\anc. 
These two sets are known as the diastolic and the systolic branches of the reduced slow manifold, 
respectively. 

To formulate the CV restitution boundary value problem (2) for equations (3), we look for 
solutions in the form of waves travelling with a constant velocity c and a fixed shape i.e. we 
assume the travelling wave ansatz z = x — ct, which gives the asymptotic boundary- value problem 

y(PE dE , , 
-TT + ^^T^ + FE{E,y)^0, (4a) 
dz^ dz 

c^ + Fy{E,y)=0, (4b) 



E{0) ^ E{P) ^ Eo, ^ 
dz 



z=0 



dE 
dJ 



, y{0)^y{P). (4c) 

Z = P 



We first formulate the slow and fast subsystems corresponding to this problem, and after that we 
will discuss matching and boundary conditions. 

The slow-time subsystem is obtained immediately from equations (4) in the limit e — > 0, and 
has the form 

FEiE,y)^0, (5a) 

c^ + Fy{E,y)=0, (5b) 
dz 

so a unique solution is obtained by imposing a a single boundary condition, e.g. 

y(0) - y*, (6) 

where is a constant. Here we use E{z) = lim^^o E{z) and y{z) = \\in^^oy{z) to denote the 
slow-subsystem solution approximation and distinguish it from the exact solution. 

The fast-time subsystem is obtained from equations (4) by first rescaling the traveling wave 
coordinate. Z ~ (z — z^)/e, where = const is the position of the jump (front or back) in the 
slow wave coordinate, and then taking the limit e — 0, which gives the equations 

d'^V dV 

^ + c-+FEiV,Y)=0. (7a) 

dY , 



the boundary conditions for which can be taken in the form 

l/(-oo)=£z, V{+^)=Er, V{Q) = Vo, (8) 

dy 

— -0, r(-cx3) = n. 

Above we have introduced V{Z) = limj^o E{z^ + eZ) and Y{Z) = Unie_>o y{zsr + tZ) for the 
fast-subsystem approximation to exphcitly distinguish it from the solution in the original slow 
coordinate. The arbitrary constant Vq is assumed in the range of V , and is used to define the 
position of the jump in terms of the slow wave coordinate z, so that E{z^) = Vq. Note that 
equations (7) are obtained in that form without the need of e-dependent scaling of the speed c 
only if the spatial scaling depending on e is chosen as in equation (3), which is the reason for that 
choice. 



3.2 Solution of the fast subsystem 

Due to equation (7b), y is a first integral, and then equation (7a) together with the boundary 
conditions (8) present an eigenvalue problem for the profile V{Z) and velocity c of a trigger wave, 
depending on F = as a parameter. It also depends, of course, on the values of the voltage to the 
left and to the right of the front, Ei and E^, which should be the two stable roots of Fe{-, Y^), i.e. 
{El, Er} = {E_{Y^), E^{Y^)}. Under the assumptions made about function Fe{-, Y^) and with an 
appropriate choice of the pinning value Vq, the fast-time boundary- value problem for the trigger 
wave has a unique solution, which is guaranteed by a result due to Aronson and Weinberger [3, 
Theorem 4.1]. We denote this unique solution by 

ViZ)^ViZ:Y,,Ei,Er); (9) 

and the corresponding propagation speed by 

c^C(Y,;Ei,Er). (10) 

The boundary value problem (7b), (8) is invariant with respect to simultaneous transformation 
Z —Z, c — c, El -H- Er- Hence, it follows that 

Viz- Y,,Er, El) = Vi-Z; Y„Ei,Er) (11) 

and 

C(Y,;Er,Ei) = -C{Y,;Ei,Er). (12) 

Note that these formal solutions can be with positive as well as negative c; we are, however, only 
interested in the waves propagating rightwards, c > 0. Now we can discuss fronts and backs as 
two different types of trigger waves. 

• Suppose that for some Yf we have C(Yf; E+(Yf), E^{Yf)) > 0. This means that we have a 
forward propagating trigger wave that switches the system from the lower quasi-equilibrium 
Er = E-{Yf) to the upper quasi-equilibrium Ei = E+(Yf). We will call this type of fast 
solution a front. 

• Now suppose that for some Yf, we have C(Yfc; E^{Yi,), E-{Yi,)) < 0. This means that an up- 
jump trigger wave does not propagate forwards but retracts backwards, and is not suitable 
for us as we are interested in forward propagating waves, c > 0. However, due to (10), we 
know that we then have C{Yb; E-(Yb), E-i-{Yb)) > 0, that is there is a forward propagating 
down-jump trigger wave switching from the upper quasi-equilibrium Er = E+(Yb) to the 
lower quasi-equilibrium Ei = E-{Yb). We will call this type of fast solution a back. 
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3.3 Solution of the slow subsystem 

Equation (5a) implies that E = E±{y) or E = E^:{y). By assumption A2, the latter solution is 
the unstable branch of the reduced slow manifold Fe{E, y) = while E±{y) are the stable ones. 
Hence ignoring the possibility of "canard" solutions that involve the unstable branch, we must 
solve 

c^^+Fy{E±iy),y) = 0, (13) 

which is separable and can be easily integrated, giving the (spatial) length of the piece of a solution 
say between y{zi) = yi and y{z2) = y2 as 

where the plus subscript refers to the systolic (action potential) branch and the minus subscript 
refers to the diastolic (diastolic interval) branch. Naturally, Z2 > z\ requires that y\ — y2 and 
Fy{E±{y) have the same sign. 

3.4 Matching 

A period of a steadily propagating periodic pulse train, to which the asymptotics described above 
are applicable, must include at least one fast front, one fast back, one systolic interval and one 
diastolic interval. The restitution curve sought for can be obtained from conditions of matching 
of these four pieces. 

The asymptotic matching of the fast and slow pieces in the leading order in e is rather straight- 
forward. Let us consider a fast jump solution {V{Z),Y{Z))^ existence of which is guaranteed by 
the Aronson- Weinberger theorem. [3, Theorem 4.1], located at z = so that Z = {z — 2:*)/e, 
and compare it with the slow solutions {E{z), y{z)) ahead and behind it. By van Dyke's matching 
rule, we have 

lim V{Z) = Ei= lim E{z), 
lim V{Z) = Er= lim E{z), 
lim y{z) = = lim J/(z), (15) 

2— >2,— 2— >Z,+0 

that is, y is continuous across z = and the jump of E at z, is related to y(z*) and c via the 
speed equation (10). 

Note that above we distinguished between y, Y and y only in order to demonstrate explicitly the 
decoupling of the slow- and the fast-time problems. However, they all coincide in the leading order 
in e, below we will, for simplicity of notation, use y to represent any of them, as this distinction 
does not run any deeper. This applies, in particular, to y, = F,, y/ = Yf and yt — Y\,. For 
the same simplicity, henceforth we write E instead of E as they coincide in the same limit almost 
everywhere. We keep distinguishing F, though, as it describes gradual change of the voltage where 
E has a jump. 

Let us take the conduction velocity c as the parameter, i.e. construct the periodic pulse prop- 
agating with a given speed c > 0, and then calculate its temporal period P. Further to assump- 
tions A1-A4, we make the additional assumption. 

A5. Function C(%j\ £'+(?/), E^{y)) is a monotonically decreasing function of y. 

This is easily verified for the two examples that follow; note that a monotonically increasing 
function can be deal with in just the same way. Then equation 

Cfiyf)=C{yf;E+iyf),E^{yf))^c (16) 
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may have at most one solution for first integral parameter yf of the front, and equation 

Ci,iyf)=C{yb;E^{yi,),E+{yf,))^c (17) 
may have at most one solution for the first integral parameter yi, of the back, and we always have 

yb>yf, (18) 

as positive values of a monotonically decreasing function are achieved at smaller values of the 
arguments than negative values. 

Now we can formalize assumption A4 in the following way. 

A4. There is a nonempty interval of c values, Xc — (cmin, Cmax), which is the interval of interest, 
such that Cb^\lc) C I'y and C/"^ (Ic) C I^. 

So, under the assumptions Al~A5, for every c G Xc there are exactly two types of fast jump 
solutions, 

• a front, at y = yf, with a pre- front voltage and post-front voltage Vj , 

Vf{Z)^ViZ;yf,V^,V^), 
Vi^E_{yf), Vj^E+iyf), yf^Cf-\c), (19) 

• and a back, at y = yb, with a pre-back voltage and post-back voltage V^, 

Vb{Z) = V{Z;yb,V^,V^), 

V^^E+iyb), V^^E^iyb), yb^Cr\c). (20) 

Let us now consider the slow pieces. Since there is only a unique choice of the y-values they 
can have at their ends, namely yb and yf as shown above, there are only two possibilities: a slow 
piece that has yb on the left (back) end and yf on the right (head) end, and vice versa. For a piece 
with yb on the left and yf on on the right, due to inequality (18) and assumption A3, equation (14) 
gives positive length of the piece only if it is a systolic piece. The temporal duration of such piece 
is the APD, and equals 

Vh 

APD = / , (21) 

J -Fy{E+{y)) ^ > 

Vf 

Similarly, for a piece with yf on the left and yb on on the right, due to inequality (18) and 
assumption A3, equation (14) gives positive length of the piece only if it is a diastolic piece. The 
temporal duration of such piece is the DI, and equals 

DI=/_^. (22) 

Vf 

Hence, we have demonstrated that in the assumptions made A1-A5, for every c G Ic there is 
exactly one, up to translations along the z axis, solution of each of the following four kinds: a 
front, a systolic slow piece, a back and a diastolic slow piece, and they can be matched only in 
a unique order. Hence for every c we have a periodic solution, each period of which consists of 
exactly one piece of each kind. 

To summarize, for every c e Xc, in the leading order in e, the temporal period of the solution 

is 

Vb 
Vf 

where j// = C/~^(c) and yb = Cb~^{c) are the unique solutions of equations (16) and (17) respec- 
tively. 
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Figure 1: (color online) CV asymptotics in the Barkley model. Parameters: a = 0.7, b = 0.1. (a) 
Schematic of the asymptotic pulse solution in the {y, E) plane. Note that in this and the next fig- 
ures, the direction of the axes {E vertical, y horizontal) is different from the traditional [y vertical, 
E horizontal), which is to comply with subsequent ionic gate figures where the transmembrane 
voltage (corresponding to E here) is on the vertical axis, (b) CV restitution curves for various 
values of e. The curve e = is the asymptotic given by (32). 



3.5 Special case: cubic fast dynamics 

Problem (7a), (8) has an explicit solution in two popular special cases, for a cubic [65] and for a 
piece- wise linear [40] dependencies i^s (•;?/). The two simple examples that follow fall in the case 
of cubic nonlinearity, 

Fe{E, y) = -A{E- E^) {E - E,) {E - E+) , 

E±,^ = E±,^{y), A^A{y)>Q, (24) 

which evidently satisfies assumption A2 as long as assumption Al holds. Assumption A5 imposes 
obvious constraints on the functions E±{y) defining this nonlinearity. 

In this case it is convenient to choose Vg = {E- + E^)/2 = {Ei +Er)/2, and then the solution 
to problem (7), (8) is 

V^Va + iEr- El) tanh (^{A/Sy/^iEr ~ Ei)Z^ /2, 

c = C{yf- EuEr) = (2A)i/2 {Vq - E,) , (25) 
and we remind that {Ei.Er} = {^^-(y/), i?+(j//)}. 



4 Asymptotic restitution curves in the Barkley model 
4.1 The model 

The functions Fe{-, ■) and Fy{-, ■) of the Barkley model [5] are given by 

FE{E,y) = E{l-E){E-^^^, (26a) 
Fy{E,y)^E-y, (26b) 

where a and h are parameters, satisfying 6 > 0, 26 < a < 1 + 6. The equation of the reduced slow 
manifold FE{E,y) = is trivial to resolve and yields the branches 

E^ =0, E, = {y + b)/a, E+ = 1. (27) 
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With this choice of the branches, assumptions Al and A2 are satisfied for y E ly = {—b, a — b). 
However, assumption A3, specifically the condition Fy{E^{y),y) < 0, narrows this down to y G 
I'y = {0,a~ b). The phase portrait and the N-shaped form of the reduced slow manifold of the 
Barkley model is illustrated in figure 1(a), with an example of a trajectory corresponding to a 
traveling wave train. 



4.2 The fast subsystem 

Substituting (27) into equations (25), (16) and (17) gives the front velocity 

c = Cfiyf) = ^(^l-2yi±^^>0, yfe{0,a/2-b), (28) 

and the back velocity 

c = Cb{yb) = ^(^2^^^-l^ >0, yte{a/2-b,a-b). (29) 

As the front and the back have the same speed c, we can obtain yb for a given y/ by eliminating 
c from system of equations (28) and (29) and resolving it with respect to yb, which gives 

yb^a-2b~yf, (30) 

and provides a link for matching with the slow-time problem. The resulting interval of achievable 
speeds is = (0, (1 - 2b/a)/V2). 



4.3 The slow subsystem and matching 

Evaluating expression (23) along the stable branches E^{y) and E-{y) of the reduced slow mani- 
fold given by (27), yields the temporal period of the wave 



P = In 



(1 - yf)yb 



(1 - yb)yf 



(31) 



Combining expressions (28), (30) and (31), finally, yields the CV restitution curve in explicit 
form 



P = In 



' (a - 25 + ac^/2) (2 - a + 26 -f acy/2) " 
{a -2b- acV2){2 - a + 26 - acV2) 



(32) 



which for c G Ic gives the range P € (0, oo). Figure 1(b) illustrates this result in comparison with 
curves obtained by numerical solution of the full boundary value problem (2) for equations (3) 
with right-hand sides given by (26) and periodic boundary conditions as described in section 2. 



5 Asymptotic restitution curves in the FitzHugh-Nagumo 
model 

5.1 The model 

We will use the right-hand sides Fe{-,-) and Fy{-,-) of the FitzHugh-Nagumo equations in the 
following form, 

FEiE,y) = Eil-E){E-p)-y, (33a) 
Fy{E,y)^aE-y, (33b) 

which is related to the original [22, 47] formulation by an affine transformation of the variables, 
involving the small parameter e. Parameters a and j3 are assumed to obey 0</3<l/2,a> 
(1 — /3)^/4. The corresponding phase portrait is illustrated in figure 2(a), with a typical trajectory 
corresponding to a traveling wave train. 
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Figure 2: (color online) CV asymptotics in the FitzHugh-Nagumo model. Parameters: /3 = 0.13, 
a = 0.37. (a) Schematic of the asymptotic pulse solution in the (y, E) plane, (b) CV restitution 
curves for various values of e. The curve e = is the explicit asymptotic result found from (35) 
and (38), (39) as detailed in the penultimate paragraph of section 5. 



5.2 The fast subsystem 

To use the general results on the front velocity (25), we need to know the branches of the reduced 
slow manifold E±{y) and E^,{y) as functions of the slow variable y. However, unlike the case of the 
Barkley model, here this would require using the formula for the roots of a generic cubic equation. 
This is rather inconvenient so we employ an alternative strategy. Given the value of the prc- 
front voltage at the lower branch of the reduced slow manifold, V/ = E_{yf), we determine from 
equation (33a) the corresponding value of slow variable during the front yf = V^{1 — V^) {V^ — fH). 
Further, to find the corresponding values of E^,{yf) and the post-front voltage of the up-jump 
= E^{yf), we need to solve the cubic E [1 — E) {E — (3) — y f =^ Q for which we already know 
one root, namely E = , so the cubic is divisible by {E — V^). Hence the other two roots are 
solutions of the resulting quadratic equation, which leads to the required values 

(y/) - ^ (/? + 1 - - ^J{P-ir + 2Vlil3 + l)-3ivIr^ , (34a) 

= l(^li + l-V^ + v/(/3-l)2 + 2Kf(/3 + l)-3(yi)2^ , (34b) 

and therefore we get the expression for the front velocity 
V2 



Cf{y) 



31/i + 3^J i/3 - 1)2 + 21// (/3 + 1) - 3(y/)2 - /3 - 1 



(35) 



assuming V/ = E^{y), and a similar expression for the back velocity. Finally, from the condition 
Cbivb) = Cf{yf) wc find the prc-back voltage as 

V^ = E+{yb)^'^{P + l)~Vl (36) 

which provides the link for matching with the slow-time system. The interval of bistability required 
by assumptions Al and A2 in this case is ly = (ymin, 2/max), where ymax,min = (1 + /? ± (t)(2 — /3 =F 
(t)(1 — 2/3 ± cr)/27 with upper signs are for ?/max, lower signs are for j/min, and a = \Jl — (3 + (3'^ . 
Assumption A3 narrows this to y € I'y = (0, j/max)- Assumptions A4 and A5 for this interval are 
verified by direct elementary calculations, giving Xc = (0, Cmax) where Cmax = (1 ^ 2/3)/-\/2. 
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5.3 The slow subsystem and matching 

To use the coordinate E± to describe the motion along the reduced slow manifold, we rewrite (5) 

as 

dz dy dz -SE±^ + 2{/3 + 1)E± - (3 

Therefore, we have the action potential duration as the time between Vj and along the upper 
branch of the reduced slow manifold, 

APD= r"^^^!±^(^±i^di., (38) 
Jyj aE-E{l-E){E-P) ^ ' 

and the diastolic interval as the time between and Vl along the lower branch of the reduced 
slow manifold, 

_ -iE^ + 2{(3 + l)E^I3 

""^-k o.E-E{l-E)[E~P)'''' 

and hence the period of the wave P = APD + DI. Note that equations (38) and (39) have the 
same integrand, and only differ in the integration limits, which are related by relationship (34b), 
a similar expression relating and Vj^, and equation (36). 

To summarize, equation (35) gives the wave velocity c as the function of the pre-front voltage 
Vf. Equation (36) gives the pre-back voltage as a function of the pre-front voltage V/. 
Equation (34b) and its analogue for the back give the post-front voltage and post-back voltage 

as functions of the pre-front voltage 1//. Using those, finally, equations (38) and (39) give 
the wave period P, as a function of the pre-front voltage V/. Hence, we have a parametric 
description of the conduction velocity restitution curve, (^'(Vj/), c(V/)) in a parametric form 
with parameter the pre-front voltage . The parametric representation can be transformed into 
explicit representation by noting that expression (35) is equivalent to a quadratic equation with 
respect to the pre-front voltage V/ and can be easily solved to give the desired explicit expression 
for P = APD -I- DI as a function of c; the result, however, is rather lengthy and we omit it here. 

Figure 2(b) presents a comparison between this explicit asymptotic P{c) dependence and the 
solution of the full periodic boundary- value problem at various values of e. 



6 Asymptotic restitution curves in the Caricature Noble 
model 

The classical asymptotic theory of slow-fast systems described in the previous sections is not 
appropriate for the asymptotic reduction of cardiac equations which have a different nature, as 
pointed out in the Introduction. To develop a fully fledged alternative general theory is beyond 
the scope of this paper. Instead, in this section we study an archetypal "caricature" model of 
cardiac excitation previously proposed in [9] . One can think of this caricature model as a simple 
example of an ionic cardiac model in which a small parameter has been embedded so as to reveal 
explicitly the non-Tikhonov properties of the equations. The resulting fast and the slow problems 
have analytical solutions in closed form which makes the model convenient for investigation of 
well-posedness of the asymptotic reduction of the CV restitution problem in this particular case. 
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6.1 The model 



We consider the following set of equations [9], 



— = -GNa (-BNa - E) e{E - E,) h + g^iE) + G{E) + e^, (40a) 
^±^\FnXB(E^-E)-h), (40b) 
^=F^{e{E-E^)-n), (40c) 

where 

92iE) = g2iO{E^ - E)+g22e{E - E^), 

921 = -2, 322 = -9, 

r h{Ei^E), Ee (-oo,£;t), 
G(£;) = <^ k2{E-E2), Ee [E^,E,), 
[ hiEs-E), Ee [E,,+oo), 

fci = 3/40, fc2 = 1/25, fc3:=:l/10, El ^-280/5, 

E2 = iki/k2 + l)E^ - Eiki/k2 = -55, 

Es = (^2/^3 + 1)^* - E2k2/k3 = 1, 

Fh = 1/2, F„ = 1/270, 
^Na = 40, -Bt^-SO, ^* = -15, Gjva = 100/3, (41) 

and where 9{-) is the Heaviside unit step function. 

The time in this model is measured in ms and the voltages E, Ei, E2, £^3, i?*, E-^ are measured 
in mV. Correspondingly, the units of g2, (721, 322 are mV/ms, and the units of Gnu, and Fn 
arc ms~^. As discussed above in section 2, the space scale is chosen to get the convenient value 
e for the coefficient at the voltage diffusion term, so the dimensionality of x in (40) is given by 
the "space unit" su = ms^/^. The real physical lengths are given by xD^/^ where D is the tissue 
voltage diffusion coefficient in the direction of wave propagation. The rest of the quantities in (40) 
are dimensionless. 

This system is obtained from the authentic Noble model of Purkinje fibers [49] using a set 
of verifiable assumptions and well defined simplifications as detailed in [9] . The main features of 
equations (40) which make them an appropriate illustration are: 

(a) They reproduce exactly the asymptotic structure of the authentic Noble model [49], which is 
guaranteed by the embedding of the artificial small parameter < e ^ 1. The authentic Noble 
model is the prototype of all contemporary voltage-gated cardiac models, and we believe that 
the asymptotic structure of (40) is rather generic in this class. Realistic voltage-gated cardiac 
models do not have explicit small parameters already present in them; or, rather, they have so 
many parameters that it is not a straightforward task which of them to use for asymptotics. 
Hence we employ a procedure of embedding artificial small parameters, as discussed e.g. in [9]. 
An example of the embedding procedure appears in section 7 below, where the Beeler-Reuter 
model [6] is discussed. 

(b) Equations (40) have the simplest possible functional form consistent with property (a). Most 
functions in the right-hand side are replaced by constants as justified in [9] which allows 
analytical solutions to be obtained in closed form. This in turn makes it possible to prove the 
well-posedness of the asymptotic boundary value problem to be formulated below. 

For brevity, wc shall call this model "Caricature Noble" . 
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6.2 The asymptotic reduction of the CV restitution boundary-value 
problem 

Model (40) contains an explicit small parameter e embedded in essentially the same way as it would 
be in a realistic model. In this section we demonstrate how this may be used for simplification of 
the Caricature Noble model or, indeed, of a more realistic ionic model. 

A slow-time subsystem which describes the plateau and the recovery stages can be obtained 
immediately from equations (40), by taking the limit e — >■ 0. At time scales much longer than e, 
the second equation implies h — >■ 9[E^ — E). Hence the first term of equation (40a) is proportional 
to 9{E — Et)9{E-^ —E) = which vanishes in the limit e despite the large factor in front of 
it^. The diffusion term d^E vanishes in the same limit and wc arc left with the slow-time system, 

^=g^{E)n^ + G{E), (42a) 

— ^F,,{0{E-E^)-n), (42b) 

where rj ~ t — x/c is the traveling wave coordinate, which wc use in this section instead of our 
standard choice of z = x — ct. A fast-time subsystem of equations (40) can be obtained by 
stretching time and space, T — t/e, X = x/e, taking the limit e — > and neglecting the equation 
for n which decouples from the rest. It is useful to distinguish explicitly the functions of the old 
from the functions of the new independent variables, say E{x,t) = V{X,T) = V{x/e,t/e) and 
h{x,t) = H{X,T) = H{x/e,t/e). It is also useful to introduce at this stage the following non- 
dimensionalization (which, as noted above, is different from other sections and specific for this 
particular model) 

V - E,. I — 
« = ^ ^, i = X^H, T = F,,T, 

g^^,C = c/^,. (43) 

In these variables, the travelling wave ansatz becomes C = t — ^/C. As a result of these transfor- 
mations we obtain the following fast-time model of the wave front, 

— = e{v^ -v)-H, (44b) 

where vi^ = {E^ — E^)/{E^g_ — E^) < 0. In a periodic wave train, a front propagates in the tail of 
the preceding wave, so slow pieces described by (42) and fast pieces described by (44) alternate, 
and there is one slow piece and one fast piece per period, as opposed to two fast pieces and two 
slow pieces in classical Barkley and FitzHugh-Nagumo models. The matching points for the van 
Dyke rule are: (a) the end of a slow piece rj = P corresponds to the beginning of the fast piece 
( — >■ — oo and (b) the end of the fast piece ( ^ oo corresponds to the beginning of the next slow 
piece ?7 = 0. This situation is summarized by the following set of boundary conditions 



V[-00)=Va, V[00)=V^, — 



= 0, 

H{-oo) = 1, w(0) = 0, (45a) 



^This is an attempt to summarize briefly the essence of the non-Tikhonov asymptotics of this and similar models. 
For a more detailed treatment, see our previous publications, e.g. [9]. 
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together with 



E{0) =E,+v^ {En, 

E{P) ^E^+Vc. {En: 

n{0) = 



a 



a 



E.), 



(45b) 



where C e (0,cxj), P G (0,oo). Va S i~~oo,v^) and G (0,1) are parameters to be found. 
Condition v{0) = is a pinning condition as discussed above in (2), i.e. we choose Eq — E^. 
Condition _ff(— oo) — 1 follows from matching condition _ff(— oo) = h{P)^ by noting that h in the 
slow time system is given by ft, = 9{E^ — E) and E{P) < Ef. The asymmetry of the conditions 
imposed at C ~^ ±oo can be understood by analysing the C — > ±oo asymptotics of the linearized 
problem: the condition on the ^-derivative, v'{—oo) = 0, and the condition iJ(oo) = are satisfied 
automatically for open sets of solutions, whereas the the condition on the (^-derivative, i''(+oo) ~ 0, 
excludes solutions exponentially growing as C ^ +oo and the condition H{—oo) = 1 excludes 
solutions exponentially growing as — —oo. 

Equations (42) and (44) together with the boundary conditions (45) form a set of coupled 
boundary value problems representing an asymptotic description of CV restitution. The slow 
system is of order 2 while the fast system is of order 3 and there are 4 unknown constants namely 
C, P, Va and v^. Hence 9 conditions are needed to select a unique solution, while (45) provide only 
8 conditions. Hence a one-parameter family of solutions may be found where the wave velocity is 
a function of the wave period, C{P). 

The asymptotic boundary- value problem (42), (44) and (45) of CV restitution is essentially 
simpler than the full one (2) for equations (40). Indeed, the small parameter has been eliminated 
and the resulting system is no longer stiff. Furthermore, the right-hand sides of equations are 
simpler and each stage of the action potential is modeled asymptotically by a system of lower 
dimension. However, to be useful the coupled asymptotic boundary value problem must satisfy 
two essential requirements: (a) the coupled problems must be well-posed (b) their asymptotic 
solution must provide a good approximation to the solution of the full non-asymptotic problem. 
It is not obvious that the asymptotic formulation of the CV restitution problem satisfies either of 
these requirements in the non-Tikhonov case under consideration. While a proof of properties (a) 
and (b) in the case of any arbitrary voltage-gated cardiac model is beyond the scope of this paper, 
in this section we prove the well-posedness of the archetypal Caricature Noble problem (42), (44) 
and (45). The convergence of the asymptotic and full solutions is demonstrated numerically. 

6.3 The fast subsystem 
6.3.1 Exact solution 

To solve the fast-time equations (44) and (45a), we follow the ideas presented in [28, 55]. Since the 
right-hand-side of equation (44a) is a piece- wise function of voltage, we distinguish three intervals 
in terms of voltage separated by and 0, or alternatively in terms of the wave coordinate ^ 
we use the intervals C £ (— oo,k], C G [k,0] and ^ G [0,cx)) with internal boundaries ( = k 
and C = for which the equations v{k) = v-^ and v(0) = are satisfied, and impose natural 
continuity conditions at the internal boundaries. Exact analytical solution of the fast system can 
be obtained by first solving the _ff-equation (44b) which is separable and independent of v and 
then substituting its solution in the voltage equation (44a). In the first two intervals, C G {—oo, k] 
and C G ['*iO], equation (44a) is then readily solvable. The internal boundary point k can be 
obtained by matching the solutions in these two intervals. To solve the voltage equation (44a) in 
the third interval G [0, oo) we use the auxiliary change of variables. 



2CV^exp((A«-C)/2), 




(46) 
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and obtain a modified Besscl equation of order C^, 



dw 
'd7 



(C^ + s^)w = 0, 



(47) 



the solutions of which are a linear superposition of the modified Bessel functions Ic^ (s) and Kq2 (s) 
of order [1]. The requirement of boundedness of the solution at infinity eliminates the K(j2{s) 
term. The value of the post-front voltage v,^ is obtained as the limit of the expression for the 
voltage as ^ ^ oo and using formula [1, (9.6.7)]. In summary, the exact analytical solution of 
equations (44) is ^ 



1, C e (-oo,^:] 

exp (k — C) I C •= oo) 



I exp ( K 

( (wt ~ «a) exp (C2(C - n)) + Va 



v{0 



1 - exp(C^C/2) 



for C e (-00,0] 
/c2(2CVgexp((Ac-C)/2) 



for C e [0, oo) 



where the internal boundary point k is given by 



5 

V, 



^ <0 



(48a) 



(48b) 



(48c) 



(49) 



and the post-front (peak) voltage is 



lim v(() 



(so/2)' 



r(C2 + l)/cJso 



(50) 



where r(-) is the Gamma function [1]. The dispersion relation is then found from the continuity 
of the derivative of voltage v a.t ( = and takes the form 



(1 - 2w„) 



so = 2C^(l- 



Vry 



1/(2C^) 



(51) 



where the prime indicates a derivative with respect to the argument s. At a given pre-front voltage 
Va , the wave velocity C of the travelling impulses can be found as a solution of the dispersion 
relation (51) and we remind that for comparison with numerical results the wave velocity should be 
transformed back to the original variables, c = C\/Fh- Once again, note that the pre-front voltage 
Va cannot be found from conditions (45a) alone and so the entire front solution is a one-parameter 
function as expected. 

The existence of solutions to the fast-time boundary value problem thus ultimately depends on 
the existence of solutions to the transcendental dispersion relation (51). Solutions are guaranteed 
by the following 

Proposition 1 For every set of parameters such that < 7J| < and C > 0, there exists a 
unique value oj the excitability parameter g = L(C,Vo,,v^) > which solves (51). 

Proof 1° We will need an important property of modified Bessel functions: the ratio 



^7+1 



(P)/A(P) 



^At V| = 0, K = 0, this solution coincides with that of Hinch [28] at n^^eff = 0- 
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is a strictly increasing function of its argument p S (0, oo) for any order 7 > 0. This follows 

directly from the estimate r'^{p) > 0, [2, p. 243]. 

2° Moreover, it is easily established from the asymptotics of Bessel functions that lim r^(p) = 

p— >o 

and lim r-.(p) = 1. 

p— J-oo 

3° The right-hand side of (51) is a composite of the function 

F^:p^pl'^{p)/I^ip), 
further depending on 7 = as a parameter, and the function 

S : so = 2CVff (1 - Wf/wa)^ 

further depending on C, Va and W| as parameters. In the assumptions made, S is obviously 
strictly increasing as a function of g, and maps (0, 00) — > (0, 00). Using the recurrence relations 
[1, (9.6.26)], we can rewrite the definition of the function as 

F-,ip) = -f + pr^{p). (52) 

By 1° and 2°, we have that F^{p) is strictly increasing (and therefore invertible) and maps (0, 00) — >• 
(7,00). Overall, we conclude that the right-hand side of (51) is a strictly increasing function of g 
defined for all 17 > and with the range of (7, 00). 

4° The left-hand side of (51) does not depend on the parameter 7 and, since the pre- front 
voltage f a < by assumption, it lies within (7, 00) which by 3° is the range of the right hand 
side. Hence a solution g > always exists. Moreover, since by 3° the right-hand side is strictly 
monotonic, the solution is unique. ■ 

Denoting by G-y the inverse function to Fj at the constant order 7, the existence of which has 
just been established in 3° above, the solution of (51) can be written as 

1 / \ -l/(2C^) 

g = L{C,v^,v^) = — (l-^j Gc2{{l-2v^)C'). (53) 



6.3.2 Bounds and asymptotics 

We have demonstrated the existence of solutions to the fast subsystem (44) and (45a). We will 
now consider some estimates related to this solution which will lead to convenient explicit approx- 
imations of the propagation speed and the minimal excitability required for wave propagation. 
We treat £ (—00, 0) as a parameter characterising the system (and omit dependence on it in 
the function notations), while C and Va as variables characterising a particular front solution. 



Lower bounds on the excitability g From 1° and 2° we know that r-y{p) < 1, where the 
inequality becomes approximate equality for large p. We shall denote this as r^{p) < 1. With 
account of (52), this implies that Fj{p) < 7 +p. A sharper upper bound, r^(p) < {{p^ + "J^Y^^ — 
7)/p, is given in [2, (9)] and implies 

F^ip)<il^+PY^^ 7>0, p>0. (54) 
Substituting this into (51), we get the more easily tractable approximation 

C^{1 - 2«„) < (4C'g (1 - v^/v^f" + C^) . (55) 
Resolving this with respect to the excitability parameter g, we get a lower bound for it in the form 

9>UC,v^)^v^{v^-l)C^^^) . (56) 
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The minimum of i(C, Va) defined in (56) with respect to the wave speed C € (0, oo) at a constant 
prc-front voltage Va G {—oo,v^) is achieved for 

C = C*M = fin f— ^ (57) 



ft 



and is equal to 



L*M=ev^{va.~l)\n[ —]. (58) 

Hence, for any fixed value of the pre-front voltage Va G (— oOjVj) and the excitability parameter 
9 > L.*iva), there exist two solutions for the wave speed C, namely Cfustiva, g) > C*(va) and 

Similarly, the minimum of L{C,-) with respect to the pre-front voltage € (— oo,t;|) for a 
fixed value of the wave speed C € (0, oo) is defined by 



Va = Va (C) = ^ — , (59) 

a^Vf~2C^Vf-C'^, 

which is the negative root of the quadratic equation 

2C^Va^ - (wt + 2C2wt)w„ +vf + C^vf = 0, (60) 

and a corresponding value L*{C) = L{C, Va*{C)) (the explicit equation for which is lengthy and 
of no further consequence so we omit it). Hence for any given wave speed C, and excitability 
parameter g > i:*(C), we have two solutions for the pre-front voltage Wq, separated by the value 

VaHC). 

Alternatively, equation (60) can be resolved with respect to C to obtain 

C = C*{va) = ( ^-^^^P^) . (61) 

^ \{1 - 2Va){Vf - Va) J ^ ' 

The absolute minimum of L(C, Va) over {(C, Va)} = (0, oo) x (— oo, v^) is achieved when C*{va) = 
C*{va) or, equivalently, 

(g:)^ ('.-';)('-^".' i„(^).i. 

\C* J -V^[l-Va) \V^-VaJ 

The left-hand side f{va) of this equation is a strictly decreasing function in Va G (— oo,W|), 
which is established e.g. by differentiation and using the estimate ln(x) > 1 — 1/.t for x > 1 
(the calculations are tedious but elementary and we omit them). Besides, lim f{va) = 2 and 

lim f{va,Vi) = 0, hence the above equation for the minimum always has a unique solution, 

Va ~ Va**, with corresponding C = C** and L = L**, all depending on Vf G (— oo,0) as a 
parameter. 

To summarize, we have established that the lower bound L{C',Va) given by (56) has a unique 
absolute minimum L** — L{C** ,Va**) in {{C,Va)} = (0,oo) x (— oo,w-|-), tends to infinity as 
C — ^ and C — > oo uniformly in Va G (— oo,W|) and as — >■ — oo and Va — > uniformly in 
C G (0, oo). Hence by [42, Theorem 3.1], all level sets L{C,Va) = const > L** are diffeomorphic 
to each other, and thus are simple closed curves circumventing (C**,Uq**). Moreover, for any 
g > L** , there are two values of Va = Wq_ and = i'q+ > i'q_, that are exactly two solutions 
of L^*{va) = g, such that for any Va G (wq„,Wq_|_), there exist exactly two solutions for the front 
velocity, C = Cfast(wa,ff) and C = Csiow{va,g), where Csiow{va,9) < C*{va) < Cfast(wa,ff)- 

Consequently, the level sets L{C,Va) = g, which are the dispersion curves defined by (51), are 
located wholly inside corresponding level sets of the lower bound L(C,Va) = g- In particular, 
there are no solutions for (51) for g < L** . 
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Estimates of the propagation velocity Inequality (55) can also be explicitly resolved with 
respect to C, giving double-sided inequality 



where 



C<C<C (62) 



1/2 

C=\ , 1"(1-V^^) \ ^ (63) 



W_i(i^^(^ln(l-t,tA'a)^ 



1/2 



Wr 



(^^-^in(i-nK)), 



and Wo(-) and W-i{-) are the principal and the alternate branch of the Lambert function [16], 
respectively. 

The upper inequality in estimate (62) becomes an asymptotic equality for the faster velocity 
in the limit of rescaled coordinate sq — >■ oo which is achieved e.g. for large excitability g and wave 
speed C at fixed pre- front voltage v^. In this limit, the arguments of the Lambert functions are 
small, and since Wa{z) = z + O {z^), z — 0, we have 

Cfast « C « , ^ , <? ^ oo (65) 

\Va[Va - 1)/ 

which agrees with [28, (39)], as should be expected since in this limit the difference between our 
model and [28] is inessential. 

Similarly, using a crude estimate for the alternate branch of Lambert function W-i{z) = 
ln(— z) (1 + (1)) for z -> —0, see [16], we find for the lower inequality in (62) 

1/2 

^ ^n^i ln(l -t't/«a) \ 

Osiow « C w — p ^ , g ^ oo. 



However, because the convergent series representation of the Lambert function W-i is in terms 
of logarithms, it converges rather slowly, and thus the asymptotic above does not give a good 
approximation for the standard parameter values. A better, more accurate and simpler asymptotic 
can be obtained directly from (51) in the limit sq — >■ and using Il{z) / I^{z) = v/ z/{2v + 2) + 
O (2:^), z — > 0, which can be easily obtained from [1, (9.6.7)]. This leads to 

Csiow^ ryr^^^^ )'^ 9->oo. (66) 



V \n{-Va/g) 

Note that since the estimate (66) is only asymptotic rather than uniform, it is not necessarily 
defined for all g for which the solution exists. Indeed, since ln(l — Wf/wa) < for all Va < v-\, 
we must require that ln(—Va/g) < in order that C S M. This is only possible when g > —Va, 
although according to (58), for V| — > 0, the minimal excitability i* — > 0. 



6.3.3 On the properties of the exact solution 

The above estimates of the propagation velocity C are derived from the lower bound of the 
excitability L(C, Va). The exact dependence L{C, Va) exhibits similar properties, as demonstrated 
numerically in in figure 3. We summarize these properties in the following 

Conjecture 1 The function L{C,Va) has an absolute minimum L'^ = L{C'^ jVa"^) in {C,Va) G 
(0,c») X (— oo,W|), and for every g > L"^ the level set L{C,Va) ~ g is a simple closed curve, 
crossing each line C = const (or alternatively each line Va = const^ at most twice. 
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Figure 3: (color online) Solutions of the dispersion relation (51). (a) Accurate numerical so- 
lutions, for the standard value oi g — 200/3 (thick solid black line) and smaller values of 
g = 18,20,30,40,50 (thin dashed red lines, g increasing from inside out) and the standard value 
of v^. (b) The accurate numerical solution for all standard parameters (solid black line) against 
the solution provided by the lower bound L{C,Va) (dashed red hue) and the asymptotics (65) for 
the fast branch and (66) for the slow branch (dash-dotted blue line). The magenta dotted lines 
indicate the position of the tip of the curve Va* as estimated by (78). 



Supposing Conjecture 1 is true, the solutions of the dispersion relation form a simple closed 
curve, for every g > L^. Then there exist and with < < v^, such that for any Va € 
(?;^,w^) equation (53) has two positive solutions for C, which we may denote C = C^{va;v-\,g) 
and C = {va;v^, g), where < C~ < C+. Since the level sets are simple closed curves the 
statement may be inverted so that C has the role of the independent parameter: for every C 
in some open interval C € (Cniin,Cmax) there exist two distinct values of the pre- front voltage 
Va which satisfy equation (53). At the end points of the interval i.e. C — Cmin and C = Cmax 
equation (53) has single solutions for v^. The points Cmin and Cmax are extremum points of C as 
a function of Va and the CV restitution curve has opposite slopes to the left and to the right of 
each of them. 

The fast-time systems is linked, via the boundary conditions (45) and parameters Va and v^j 
to the slow-time system the solution of which is discussed below. 



6.4 The slow subsystem 

The slow-time problem (42) and (45) can be solved exactly, too. Notice that it does not depend 
on the wave velocity and it is therefore similar to the slow problem considered in [9]. First, we 
solve equation (42b), which is separable and independent of E. Its right-hand side is different in 
each of the two intervals rj G [0,7?^] r] G [??tj J^] due to the presence of a Heaviside function. The 
equation is constrained by the continuity condition at 77 = yy^ and periodic boundary conditions 
at 77 = and 77 = P. The solution has the form 

77,(77) = n{r]) = + exp(-F„77), (67) 

where the values of the constants 6i and mi are different for the intervals 77 G [0, 77^] (i = 1, 2) and 
77 S [771, P] (i = 3) and arc given in Table 1, and the overset index here above 77 and below above 
E designates the corresponding interval. This solution is then substituted in (42a) which becomes 

E' + a,E = ft E (;) ^^i"''^ exp(-P„77))' + 7,, (68) 
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Table 1: Values of the constants used in the expression (67)-(69). 



where the constants ai , f3i and ji also depend on the interval and are given in Table 1 . The values 
in the table are obtained by a straightforward manipulation of the model definitions of G{E) and 
g2{E) in (41) and the binomial theorem is used for the term n{r]Y. Within each of the intervals 
[O,??*], [^*)^t] ['7t'-P]' equation (68) is a first order linear ODE with constant coefficients, and 
the solution of (42a) can be written in the explicit form, 

E{r^) = m ^^ + 9. eM-a.v) + ft E (^) m[ ^J^^tM^^ (gg) 

where 6i are integration constants. The exact solutions (67) and (69) contain seven parameters, 
namely 77*, rj-^, P, 61, 02, 03, Va, that need to be found from the boundary conditions and from 
the internal matching conditions, 

E{Q)=v^{E^,-E,)+E„ E{r]^) = E,, 

EM^E,, E{vt)^E^, (70) 

l(7?t) = E^, E{P) = (^Na -E,)+E,. 

The existence of solution to the slow-time boundary value problem ultimately depends on the 
existence of solutions of the transcendental equations (70). Based on the uniqueness and existence 
of solutions to an initial-value problem, it is obvious that the problem reduces to four essential 
unknowns, which we denote E^ = (i?Na — E^) + E^ = E{0), Ea = Va (-^Na — i?*) -I- -E* = E{P), 
n, = n(0) = n{P) and P. 

Proposition 2 Suppose that the values of the parameters of the slow subsystem (4-2) and (4-5) 
obey the same qualitative relationships as the default values, i.e. Ei < E^ < E2 < E^ < E3, 
ki, k2, k^, Fn > , 921,922 < and G{E) is continuous. Then for any e (0,1) and E^^ > E-^, 
the system of equations (67), (69), and (10) has a unique solution for Ea and P. For a fixed n*, 
this defines a function t-j- Va with domain e (wj,+oo) and range G (— oo,Wj), which is 
monotonically decreasing. 

Proof is evident from the phase portrait shown in figure 4. Every trajectory starting above 
E = Elf goes to the right and therefore eventually goes down below E-^ . Whilst below E'-f it goes to 
the left so 71(77) —>■ as 77 CO (this is also evident from the analytical solution). Therefore there 
exists a point 77 such that 71(77) = 7i(0) = tt.*. Moreover, such rj is unique, as the domain E < E-^ is 
absorbing and 71(77) is monotonically increasing outside it and monotonically decreasing inside it. 
So we have the proposed mapping (77*, Wt^) 1— >■ {Ea,P). The monotonicity of this mapping follows 
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Figure 4: (color online) The phase portrait of the slow subsystem (42) of the Caricature Noble 
model. Standard parameter values (41) arc used except F„ which increased three times to F„ — 
1/90 ms~^, for visualization purposes. Red dash-dotted lines represent vertical isoclines dn/dt = 0. 
Blue dotted lines represent horizontal isoclines dE/dt = 0. Thin solid black lines with attached 
arrows represent trajectories. 

from the fact that the trajectories cannot intersect, so if Ei^2 > ^wi, then the contour made by 
the straight line between points (n^jEai) and {n^,,Eaji) and the segment of trajectory joining 
these two points, lies within the similar contour made by points (ri*,Eai) and {n*,Euji). m 

6.5 Matching and well-posedness 

The CV restitution curve can be obtained by combining the results of the slow and the fast 
subsystems. According to Proposition 1 and Conjecture 1, for sufficiently large values of the 
excitability g the fast subsystem defines the wave velocity as a function of prefront voltage, C = 
C^{va), for Va G {—oo,v-^), which, via equation (50), defines the post-front voltage = v^^{va). 
On the other hand, according to Proposition 2, the wave period P and the peak voltage E^ = 

+ Waj(-E'Na ^ £-*) are functions of E^ = -E* + VaiE^a, ~ £■*) and Hence, the matching of the 
fast and slow solutions, for any given n*, can be obtained by solving the simultaneous system of 
equations for Va and v^), one resulting from the fast subsystem and one resulting from the slow 
subsystem, the latter depending on n« as a parameter. This subsequently provides where C{va) 
is given by the fast subsystem, and P = P(ri*,Uij) from the slow subsystem. Hence, we have 
ultimately the CV restitution curve (c(n*), P(nH.)) in parametric form and we have proven that 
the asymptotic CV- restitution problem (42), (44) and (45) is indeed well-posed. 

The equations of the aforementioned system for Va and are complicated and finding ana- 
lytical solutions seems seems unfeasible. Some qualitative insight can be obtained from numerical 
analysis, which for the standard parameter values is illustrated in figure 5. The solutions cor- 
respond to the intersection of the closed fast-subsystem contour (dashed blue line) with a slow 
subsystem line (solid red lines). The family of slow-subsystem lines stretches continuously, but 
non-monotonically, from the vertical line Va = wi at = to the vertical line Va = = 1. 

In accordance with Proposition 2, these lines are monotonically decreasing except for the above 
mentioned vertical lines for extreme values of n*. Almost all of these lines intersect the fast- 
subsystem contour, with the exception of the lines with very close to unity. The line Va = wi, 
n* = corresponds to the limit P — > (X) of the restitution curve, as that limit corresponds to a 
solitary wave propagating through the resting state, which is characterized hy E = Ei and n = 0. 
The other extremity corresponds to the line which is tangent to the slow-subsystem contour 
for a value of very close to but smaller than unity, and at a value of the pre-front voltage Va 
very close to but smaller than that of the parameter v^. As evident from the above analysis, e.g. 
see figure 3(a), we have Va* — > — as g — > oo, which motivates consideration of asymptotics 
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Figure 5: (color online) Matching solutions of the fast subsystem eqcaric-fast and slow subsys- 
tem (42) of the Caricature Noble model in terms of the dimensionless pre- and post-front volt- 
ages Va and v^. Standard parameter values (41) except F„ = 1/90 ms~^ as in figure 4. The 
non-dimcnsionalized resting potential vi = {Ei — E^)/{E-^g, — -E*) ~ —1.4242 and the non- 
dimensionalized h- and n-gate switch potential w-f = {Eif — i5H.)/(iJNa — i5*) ~ —1.1818. The 
choice of values of in (a) is: 0, 0.5 and 1.0; for in (b): from to 1 with step 0.1; in (c): from 
to 0.9 with step 0.1, then to 0.99 with step 0.01, and then to 0.999 with step 0.001. 



related to this limit. The details are presented in the next subsection. 

Another qualitative feature evident from figure 5 is that the slow lines are nearly vertical for 
larger v^. This, and the non-monotonic behaviour of the horizontal (wq,) position of the slow- 
subsystem lines on around larger values of Vi^ is a direct consequence of the non-monotonic 
behaviour of the trajectories which can be observed in figure 4: a typical "tail" of a trajectory 
starting at a large is a curve which starts at wi at n = 0, then decreases and then increases up 
until vf, besides, all trajectories starting from larger join together very closely. As follows from 
the analysis in [9] , this is due to an extra small parameter present in the slow subsystem, namely 

Fn 0. 

These observations motivate consideration of further asymptotics to the obtained solution, 
which lead to less accurate but more explicit results. We present them less formally than the main 
limit e — > as they are of a secondary importance to our main results. 

6.6 Further asymptotics 

The fast branch of the restitution curve For values of 1 — n ~ 1, the period P is large 
compared to the parameter Fn^^ which plays the role of a time constant in the n-gate equation 
(40c) , and hence we can exploit the Tikhonov singular perturbation in terms of Fn understood as 
a small parameter. This corresponds to the secondary asymptotic embedding 62 — > considered 
in [9]. In this limit, the trajectories differ only in the initial post-overshot stage, after which they 

/ ~ \l/4 

move along the reduced superslow manifold n « ^{E) = I —G{E)/g2{E) J with the exception 

of the repolarization stage when M{E) ^ M. It is important to note that except during the the 

initial transient, the trajectories are nearly the same, up to a correction which is exponentially 

small in Fn as evident from figure 4. We consider two parts of a typical trajectory, one for 

77 G [0, 77t] when n increases, and the other for r] G [j^t?^] when it decreases. Due to the above 

mentioned convergence of trajectories, the value rimax = '^('?t) practically independent on E^ up 

to exponentially small corrections. A typical value of ni„ax can be found e.g. in the following way: 

1 

first, consider solution (69) for i = 1 and 6*1 = and solve the matching condition E{ri^,) = E^ 

2 2 
for rjf, next determine 02 from the initial condition E{t]^:) = E^, then solve E{ri^) = E^ for rj-^, 

and finally with the knowledge of and n», the value of rimax can be obtained from (67) as 

TT-max = = = 1 — (1 — n*)e^^"''t. This however leads to a transcendental equation. 

Numerical value for the standard parameter values is rimax ~ 0.7827. 
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Using (67), the duration of the second half of the trajectory, between ~ rj^, n ~ ^maxj E = 
and r] ^ P, n ^ n^, E = Ea, is given by 

% = -P - 7?t = 4~ ("-max + (1 - '^max) 6^"^) (71) 

The evolution of E during the second half of a typical slow trajectory, is described by the relevant 
form of equation (42a) 

^ = .g2in„.ax' e-4^"("-"t) + k,{E^ ~ E) 

wherefrom 

E„ = E^ + f^"7"" e-^^-^^- + (e^-Ei- f^"7"" ) 6-^1"^ (72) 
ki — AFn \ ki — 4F„ J 

Combining (71), (72) and (65), we get an explicit dependence of c(P) in elementary functions. 

The slov^r branch of the restitution curve The slow branch is considered in a similar way. 
The difference is in the initial transient where a typical trajectory approaches the reduced super- 
slow manifold from lower values of E rather than from the higher E as it was for the faster branch, 
and also in the c{va) dependence. Hence the dependence of c(P) is obtained by combining (71), 
(72) and (66). 

The turning point of the restitution curve The turning point is the point where the fast 
branch meets the slow branch. It is characterized by extreme proximity of Va to and to 1. 
The front parameters can be estimated via the limit ^ U|, g — > oo of (57) and (58), which 
gives the highest pre-front voltage as 

v^*^v^{l + {-v^{l-v^)/gr) (73) 

and the corresponding slowest stable front velocity as 

1/2 



C* = CW).^eln^— . (74) 
Given Va* and C* , equation (50) then gives the value of the post-front voltage. 



soV2 

1 ^ (75) 

r((C*)2 + l)/(c.).(so*) 

so* = 2C* V5 (1 - ^^tK*)'^^'^'''^'^ • (76) 

The duration of the slow trajectory corresponding to these v^* and Vuj* will be an estimate of the 
shortest wave period possible in this model, P*. A simple approximation of it can be obtained 
from the consideration that in the limit Va Vj^, we have n(r]) « 1 throughout, hence the slow- 
subsystem equation for E is simplified by replacing 71(77) = 1 so the link between n and E dynamics 
is only via the values of 771 and P. The period P is almost the same as the rjf taken by the voltage 
to decrease from Ei^ to E^ , since the interval P — -q^ needed to decrease from E-^ to Ea is small 
compared to that. In this case, the interval t/i can be estimated from solutions (69) for 7 = 1,2. 
Hence, we have 

p*^l, ( -g22 + k^{E2-E^) \ 1 ( k^{E^* ^ E^) - g22 \ , - 

^ k2 \^922 + k2{E2-E,)J ^ ks \ ksiE, ~ Es) - 922 J ' ^ ' 

Equations (73)-(77), together with the scaling relationships (43) define the turning point of the 
restitution curve. For the standard parameter values, this gives 

C* w 2.973, 7>„* w -1.18199, ~ ^ 1.709 • 10"^, 

c* w 2.102 su/ms, P* w 13.09 ms. (78) 
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Figure 6: (color online) Restitution curves of the Caricature Noble model, (a) in Cartesian and (b) 
logarithmic coordinates. Insets in panel (a) show selected features magnified. Lines in all plots as 
described by the legend in panel (b), specifically, accurate numerical solution of the full boundary- 
value problem (2) for various values of e > (bold solid black, short-dashed green and dash-dotted 
cyan for e = 1,0.7,0.3 respectively), the accurate matching of the fast and slow subsystems (42)- 
(45) which corresponds to the limit of e = (long-dashed red), and the approximations of the fast 
and slow branches of the RC at e = given by asymptotics (71), (72), (65) and (66) (dash-dotted 
blue). The magenta dotted lines indicate the position of the tip of the curve Va* as estimated by 
(78). 



6.7 Comparison of the asymptotics with the exact solution 

With the aim to assess the proximity of the analytical solutions obtained in the main asymptotic 
limit e — > and the full numerical solutions of the CV restitution boundary value problem, we 
present in figure 6 sets of CV restitution curves of the Caricature Noble model, and we show in 
figure 7 the action potential profiles for two selected base cycle lengths P. We also demonstrate in 
figure 6 the asymptotic estimates of the upper and lower branches and of the tip position P* of the 
curve found with the help of the secondary embedding Fn 0. The asymptotic CV restitution 
curve was obtained by solving numerically problem (42), (44) and (45) which defines c{va) and 
P{va). The full CV restitution curve was obtained by solving the full boundary- value problem 
(2) formulated for equations (40). Figure 6(a) presents the curves in Cartesian coordinates and 
figure 6(b) in semi- logarithmic coordinates to reveal in more details the behaviour at small values 
of the wave period P. We can see that as e is decreased the solution of the full problem converges 
to the solution of the asymptotically reduced problem at e = 0, and that the full model curve 
for e = 1 at standard parameter values is close to the asymptotic limit everywhere except at the 
smallest values of the period P. 

This indicates that at small P, the parameter e is not a "good small parameter". Note that 
in our asymptotic analysis, we have calculated the period P as the length of the slow subsystem 
solution, and we have neglected the contribution of the fast subsystem, i.e. the duration of the 
front, which is small of the order O (e). However, at the smallest values of P, the the front length is 
comparable to duration of the solution of the slow subsystem. This can be seen already in figure 6 
and figure 7, and is further confirmed by the analysis of the dependence of the minimal wave period 
P on e shown in figure 8. We see that to second order, the basic cycle length can be approximated 
by P = Po + ^Pi + O (e^) , where Pq is the cycle length given by the asymptotic theory presented 
above, and Pi « 13.4 may be interpreted as the front duration in the fast subsystem. Note that 
Pq ~ Pi, hence neglecting O (e) for smaller P produces relatively large error. That this is not the 
whole story, however, as not only the horizontal position of the P* = Pmin point changes with e, 
but also its vertical position c* , so a proper next-order asymptotic should take into account of the 
influence of the slow subsystem on the front velocity as well. 
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Figure 7: (color online) Profiles of the action potentials corresponding to the faster branch of the 
solution of the Caricature Noble model. The values of the period P and the values of e are given 
in the plot. The values e 7^ correspond to the full CV restitution boundary value problem (2) 
for (40) with standard parameter values while e = corresponds to the asymptotic limit given by 
(42), (44) (45). 



7 Asymptotic restitution curves in the Beeler-Reuter model 

In this section we apply the methodology presented above to the Beeler-Reuter ventricular model 
[6] . This model is an example of a realistic voltage-gated model which represents an intermediate 
step between relatively simple early models and complicated contemporary models. It has played 
an important role for understanding of cardiac excitability with a large volume of literature devoted 
to it, and it is still the model of choice in many situations, e.g. [13, 20, 27] are some recent 
examples. In the following, we find the CV restitution curve of the Beeler-Reuter model using 
both the asymptotic formulation and the full formulation of the periodic boundary value problem 
as described above and demonstrate an excellent quantitative agreement. 



7.1 The model and the asymptotic embedding 

The initial step of our approach requires an appropriate asymptotic embedding of the original 
Beeler-Reuter model [6]. The embedding is constructed following the procedures presented in 
detail in our earlier works [9, 11, 55] and here we shall summarize briefly the relevant arguments. 
We would like to remark that an analogous embedding procedure applies to the Caricature Noble 
model of section 6 where e appeared seemingly without much justification. 
We rewrite the Beeler-Reuter model in the one-parameter form, 

dE 1 0"^ E 

-g^ = -ffNa (-^Na - E) mo^(E; e)i h + I^{E, y) + e— , (79a) 

^^^\{KAE-e)-h)/Tn{E), (79b) 

% = {3UE)-j)/r,{E) (79c) 

^ = Fy(i?,y,...), (79d) 

where only the equations affected by the artificial small parameter e ^ 1 are shown. As in the 
previous model, the voltage E is measured in mV, time in ms, the gating variables h, j, m, y are 
non-dimensional, and the space coordinate x is measured in su = ms^/^. 
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Figure 8: (color online) Dependence of the minimal basic cycle length Pmin on the embedding 
parameter e. The asterisks represent Pmin calculated at the specified values of e, the solid circle 
is the estimate (78). The dashed line is the best cubical fit of all the asterisk points, which is 
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13.34 + 13.40e - 2.91e2 + 0.856^ [ms] 



The functions Ty{E) and yoo{E) are time-scaling functions and quasi-stationary values of the 
gating variables, respectively. Functions moo{E;e) and hao{E;e) are "embedded", i.e. they are e- 
dependent versions of nioo{E) and hoo{E) such that moo{E; 1) = moo{E) and h^{E; 1) = hoo{E) 
on one hand and moo {E; 0) = 9{E — Em) and hoo {E; 0) = 0{Eh — E) on the other hand, with 
Em = -33.75 mV mid = -71.33 mV so that m^{Em) ^ 1/2 and hoo{Eh) = 1/2. The last two 
parameters are analogous to E^, and E-^ of the Caricature Noble model. The rest of the model (79) 
is the same as defined in [6], namely /e = — (/Na^ + -^Ki + Ixi + Is) is the sum of all slow currents, 
y = {xi, d, /, [Ca]) is the vector of all slow variables in addition to gate j which is also slow, and 
Fy(-) stands for the functions governing the dynamics of the gating variables y. The rationale for 
this parameterisation is the following. 

1. The dynamic variable m is a 'superfast' variable and has been adiabatically eliminated by re- 
placing it with its quasi-stationary value moo- The variables E and h are 'fast variables', i.e. 
they change significantly during the upstroke of a typical AP potential, unlike all other vari- 
ables which change only slightly during that period. The relative speed of the dynamic vari- 
ables is estimated by comparing the magnitude of their corresponding 'time-scaling functions' 
as illustrated in figure 9(a). For a system of differential equations dwi/dt = Fi{wi, . . .w^), 
I = I, . . . N the time scale functions are defined as r; {wi,...) = \dFi/dwi\ \l^l...Nand 
coincide with the functions Ty(-) already present in (79). 

2. The dynamic variable E is fast due only to one of the terms in the right-hand side, the large 
sodium current (/Na (^-Na — E) moo{E] e) j h, whereas other currents are not that large and 
so do not have the large coefficient e'-^ in front of them. 

3. The fast sodium current gNa (^'Na — E) rn^{E] e) j h is large only during the upstroke of the 
AP, and not that large otherwise. This is due to the fact that either gate m or gate h or both 
are almost closed outside the upstroke since their quasistationary values moo {E) and hoo (E) 
are small there as seen in figure 9(b). Thus in the limit e — 0, functions iriooiE) and hoo{E) 
have to be considered zero in certain overlapping intervals E S (— oo, Em] and E € [Eh, +oo), 
and Eh < E,n, hence the representations rn^{E; 0) = 9{E — Em) and hoo{E; 0) = d{Eh — E). 

A more detailed discussion of the parameterisation given by (79), as well as the justification of our 
method of parametric embedding, i.e. a seemingly "arbitrary" introduction of an artificial small 
parameter e, can be found in [9]. 
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Figure 9: (color online) (a) On the left ordinate: Time sealing functions Ty of the dynamical 
variables of the Beeler-Reuter model during a typical action potential. The time scaling function 
of the voltage te is given by a solid line (red online), the time scaling functions of the slow variables 
are given by dashed lines (green online) and of fast variables by dotted lines (blue online) with 
all lines labeled correspondingly in the plot. On the right ordinate: Voltage of a typical action 
potential given by a dash-dotted line, (b) The quasistationary values to^ and hoc and their 
approximations rn^{E, 0) and h^{E, 0). 



7.2 The asymptotic reduction 

We are now ready to formulate the asymptotic CV restitution problem as a set of coupled boundary 
value problems similar to those described in section 6.2. Being interested in propagation with 
constant velocity and fixed shape, we introduce the travelling wave coordinates z = x — ct for 
the slow subsystem and Z = X — cT ~ z/e for the fast subsystem. As before, we distinguish the 
functions of the fast time by name and set E{z) = V{Z) and h{z) = H{Z). 

The slow-time subsystem is obtained in the limit e — >■ of the original slow independent variable 

z, 

c—+I^{E,y)^0, (80a) 
az 

^d, ^ ^ 0, (80b) 

az Tj(E) 
dy 

c^+Fy{E,y,...) = 0, (80c) 

E{z^O) ^ Ea, E{z^cP)^E^, 

3{z^Q)^:i{z^cP)^j^, (80d) 
y(^ = 0)=y(z = cP). 

The fast-time subsystem is obtained in the limit e — >■ of the fast independent variable Z , 

d^V ^ ^ .gNa(SNa - V) J 0(F - E,n) H = 0, (81a) 



dZ2 dZ 

+ (81b) 

dZ Th\V) 

V{-^) = Vo,, V'{+^) = 0, V{0)^Eh, 

V{+oo)^V^, H{~oo)^l, (81c) 

The boundary conditions of the fast subsystem include the pinning condition eliminating the 
translational invariance along the Z axis. This problem depends on four parameters, namely the 
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Figure 10: (color online) Solutions of the fast subsystem (81) of the Beeler-Reuter mode for a 
selection of values of J. (a) Front speed vs prefront voltage, (b) Postfront voltage vs prefront 
voltage. In both panels, the prefront voltage is shown in logarithmic scale with respect to Eh, 
since the curves are very close to the line Va = Eh- Note that Va increases from left to right. 



pre- front voltage Va, the post-front voltage VL,, the fixed value of the j gate inherited from the 
slow system J and the wave speed c which are determined by matching with the slow subsystem, 
given by the conditions, 

Ea=Va, E^=V^, ja = J. (82) 

7.3 The fast subsystem 

The fast subsystem (81) describes the wave front of an action potential as it propagates by dif- 
fusion. This problem has, on one hand, differential equations of cumulative order three and four 
parameters, and on the other hand, five constraints, thus the solution will typically depend on two 
"leading" parameters chosen arbitrarily, and the other parameters will be functions of these two. 
The structure of the equations (81) is very similar to the fast subsystem of the Caricature Noble 
model (44) and (45a). In fact, if tji{V) was a piecewisc constant function with a step at V = Eh, 
then the Beeler-Reuter fast subsystem would be equivalent to the Caricature Noble fast subsystem 
up to parameter values and identification of gisiaJ h^ the former with ^Na in the latter. Hence we 
may expect that the set of solutions here has a structure similar to that of the Caricature Noble 
model. In particular, we expect that J and can be determined as univalued functions of Va 
and c. Further, we expect that for a fixed J, we have the set of solutions in the {Va, c) such that 
there exists jmin such that if J > j^nin then there exists an interval {VaininiJ)i^ama.:>!-iJ)) such 
that for any Va within this interval, there are two solutions for the front velocity c = c^{J,Va), 
and correspondingly two values of the postfront voltage V^i = VL,^(J, Vq). This is confirmed by 
numerical analysis of this problem, see figure 10. 

7.4 The slow subsystem 

As in the Caricature Noble model, the slow subsystem in the leading order does not depend on 
diffusion, and therefore coincides, up to the scaling of the independent variable, with the slow 
subsystem of the single-cell model. The slow subsystem depends on four parameters, namely the 
pre- front voltage Ea, the post- front voltage E^^, the initial value of the j-gate ja and the wave 
period P, and has differential equations of cumulative order of dim(y) -I- 2. On the other hand it 
has dim(y) -|- 4 constraints, hence, similarly to the fast system, it has typically a two-parametric 
family of solutions. From the viewpoint of matching with the fast-time problem and in analogy 
with the Caricature Noble case, one possible convenient choice of leading parameters is Ea and 
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ja from which E^, P as wcU as y(0) can be found. However, unhkc the Caricature Noble case, 
it is now more difficuh to establish rigorous conditions for existence of solutions. This, however, 
can be easily done numerically. 

7.5 Matching and comparison with the exact solution 

The three constraints (82) offset the four free parameters of the slow and fast subsystems so that 
the resulting set of solutions is typically one-parametric, i.e. it is curve in the parameter space. 
The projection of this curve on the (c, P) is the CV restitution curve. Given appropriate analytical 
approximation of the relevant dependencies, which could be obtained for instance by asymptotic 
means or by fitting, solving the resulting finite (transcendental) system will produce analytical 
approximation for the restitution curve. Doing so for Beeler-Reuter model is however beyond the 
scope of this paper and we restrict to demonstrate the validity of our asymptotic approach for this 
model by solving the asymptotic matching problem numerically. 

In solving the problem (80), (81) and (82) numerically, the following features need to be taken 
into account 

(a) The fast-time problem is posed on an infinite interval. 

(b) At the same time the slow-time system is posed on a finite interval. 

(c) The length of that finite interval is the wave length of the periodic travelling wave, i.e. it is 
an unknown variable. 

(d) The fast-time system has piece-wise right-hand sides. 

(e) The pinning condition needs to be imposed in the fast subsystem. Since the fast system is 
piece-wise, it is convenient to impose it at a boundary between pieces. 

(f) The slow variable j appears as a parameter in the fast-time system. 

These features can be addressed by a number of well-known techniques and we refer the inter- 
ested reader to [4] for a general discussion and to [55] for a numerical implementation of a similar 
problem. In short, the issue of the boundary conditions at infinity can be resolved by considering 
a finite interval with boundary conditions obtained as a solution of the problem linearised about 
the asymptotic equilibria. This finite interval is then dissected into three subintervals to take care 
of the piece-wise definition of the equations. The three subinterlals together with the interval on 
which the slow-time system is posed are then mapped to the interval [0, 1] by introducing ap- 
propriate scaling factors. The pinning condition can be easily incorporated at one of the internal 
matching points. Finally, in this representation equations (81)-(82) can be solved by any standard 
boundary value problem solver such as Maple's dsolve [62], NAG's d02raf [46] and others. 

The resulting asymptotic CV restitution curve is shown in figure 11 by a dashed thin red line. 
The bold solid black line in the figure corresponds to the CV restitution curve found from the 
full non-asymptotic boundary value problem (2) written for the full Beeler-Reuter model (79) at 
e = 1. The two curves demonstrate a good quantitative agreement. 

We would like to emphasize here that, while it was still possible to solve the full problem 
numerically in this case, and the asymptotic problem was solved numerically too, solution of the 
full problem was a substantially more difficult task than the computation of the asymptotic CV 
restitution curve. The problem is certainly well-posed but it is very stiff and it required a prolonged 
experimentation with a variety of software tools and parameter continuation techniques. It is also 
worth recalling that the Beeler-Reuter model is not as complicated as contemporary models are, 
which leads us to expect that the non-asymptotic problem for such models is even more difficult 
to solve. 
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Figure 11: (color online) (a) The numerical CV restitution curves of the full non- asymptotic 
problem (2) for the original Beeler-Reuter model [6] (solid line, red online) compared to the 
numerical CV restitution curve of the asymptotic problem (80), (81) and (82) (dashed line, black 
online), (b) The AP profiles corresponding to P = 79ms and P = 500ms (same line types). 

8 Discussion 

Summary We have demonstrated that singular perturbation theory based on the largeness of 
the maximal value of the sodium current /Na compared to other currents and quality of the /Na 
ionic gates (smallness of moo{E) and hao{E) in certain voltage ranges), is capable of reproducing 
essential spatiotemporal phenomena, using conduction velocity restitution curve as the simplest 
nontrivial example involving both the fast scale and slow scale. We have explicitly compared the 
mathematical technique involved here, with similar problems in the classical FitzHugh-Nagumo 
(FHN) like models of excitable media. Apart from the different number of equations and the more 
complicated right-hand sides, we have identified in the cardiac models qualitatively new features 
of topological nature. 

Classical simplified excitable models vs ionic cardiac models Figure 12 illustrates the 
fundamental difference between FHN-like and cardiac models as far as the problem of and E^^ 
selection is concerned. In FHN-like systems (panels (a) and (b)), if the values of the slow variables 
are given, then and E^ are uniquely defined, up to the choice between front (jump up) and back 
(jump down). So, for one slow variable y as in the examples considered in section 3, there are one- 
dimensional manifolds representing all possible fronts and backs. In contrast in cardiac equations 
(panels (c) and (d), we have two-dimensional manifolds in place of one-dimensional manifolds. 
This is related to the fact that in the parametric embedding we use, the transmembrane voltage, 
like a double-faced Janus, is both a "slow" and a "fast" variable. In its slow capacity, it contributes 
to the prc-front conditions, among other slow variables. In its fast capacity, it participates in the 
excitation front. So given one "truly slow" variable (n* in panel (c) and J in panel (d)), we 
have not a one-dimensional, but a two-dimensional manifold representing possible fronts, as even 
when the value of that slow variable is fixed, the prefront voltage Ea and, consequently, post-front 
voltage E^ are still undefined. Hence to obtain the restitution curve, which is a 1-dimensional 
manifold, we have to find an intersection of the manifold representing possible fast fronts with 
the manifold representing possible slow solutions, which is something that we don't do for the 
FHN-like systems. 

Another qualitative difference is, of course, the number of fast solutions per one excitation 
pulse. In FHN-type systems, it is essential that there is a back corresponding to every front, as 
the systolic and diastolic pieces of the reduced slow manifold are separated from each other. In our 
asymptotic embedding of cardiac models, the systolic and diastolic pieces of the slow set (which is 
now not a manifold) are connected via a piece where both n and m gates are firmly closed and E 
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Figure 12: (color online) Prc/post-front voltage selection in the four different models, (a) Barkley, 
(b) FitzHugh-Nagumo, (c) Caricature Noble, (d) Beelcr- Renter. (a,b) Relationship between Va, 
Voj and corresponding value of the slow variable y in bold lines, with their projections on the 
coordinate planes in thin lines. Each of the cases has two dependencies, for the front of a pulse 
and for the back of a pulse. (c,d) Relationship between Eq,, and an indicated slow variable, as 
semi-transparent surfaces, as following from the relevant fast and slow subsystems. The lines of 
intersection of the surfaces and their projection onto {Ea, E^^) plane are also shown. 



is in its slow- variable capacity. Here it should be noted that existence of a back is not a necessary 
feature of a FHN-type system, as in presence of two or more slow variables it is possible to have a 
slow manifold with a cusp singularity so that its systolic and diastolic pieces of the are connected 
via a monostable pieces, as in the example suggested by Zeeman [64]. However, although the 
structure of ionic models admits, in principle, such manifolds [56], it has not been identified in 
any cardiac models so far. 

The numerical method for dynamic CV restitution curves In this work, we have proposed 
a computationally efficient method of calculating an ideal case of the so called "dynamic" or 
"steady-state" restitution protocol, with exactly periodic propagating pulses. The method exploits 
asymptotic splitting the problem into two parts, the slow and the fast subsystems. The advantage 
of such a split is that each of the subsystems no longer depends on the small parameter due to the 
largeness of /Na and quality of its ionic gates, so they are significantly less stiff than the original 
full problem, hence the computational efficiency. The well-posedness of the problems arising from 
the asymptotic splitting is not obvious a priori, since the asymptotic embedding we use is non- 
Tikhonov, and the general results from singular perturbation theory are not applicable to our case. 
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Wc have thus taken care to prove the well-posedness, at least for the case considered. 

Perspectives: more complicated spatiotemporal regimes There are several other proto- 
cols used in defining restitution curves, see e.g. [53] and references therein. They do not correspond 
to periodic steadily propagating wave solutions. Nevertheless, we expect that the proposed asymp- 
totic splitting should still work there, under some natural assumptions. 

So, if propagation of waves is not too complicated, then evolution of the system away from the 
fronts can be expected to be well approximated by the corresponding slow subsystem, which is 
a non-stiff ordinary differential equation for every point of the medium. After application of the 
propagating wave ansatz, the slow subsystems of the Caricature Noble and Beeler-Reuter become 
systems (42) and (80) respectively. Here we note, that without such ansatz, the slow subsystems 
would be differential equations with time as the independent variable, with exactly the same 
structure as (42) and (80) only with a different scaling, so all the above analysis applies. Fronts 
themselves for most part of their evolution can be considered locally as steadily propagating nearly 
plane waves, well described by systems like (44) and (81). 

Under these assumptions, according to the analysis of the fast subsystem, the conduction 
velocity and post-front voltage can be found as functions of the prc-front voltage and, if relevant, 
pre- front value of gate j, and the dynamics of the front is governed by an ordinary (in the case of one 
spatial dimension) differential equation for the position of the front as a function of time. This gives 
an unusually coupled system of ordinary differential equations: the local dynamics provide right- 
hand sides for the equation of motion of the next front, and trajectory of the front provides initial 
condition for the subsequent piece of slow dynamics. This unusual ODE system can predict non- 
stationary evolution of the excitation patterns, including restitution protocols. In different settings, 
this approach has been used e.g. in [10] and [21]. As demonstrated in [21], such unusually coupled 
ODE system can form a basis for investigation of such complicated spatiotemporal phenomenon as 
the spatiotemporal dynamics of cardiac alternans. ^ So extension of the present results to other, 
more complicated and important spatiotemporal regimes, seems to be a natural and imminent 
next step. 

Perspective: the problem of restitution memory The above considerations lead us to the 
problem of rate-dependence of restitution curves, or the so called "memory" effects, see e.g. [53] 
and references therein. A typical approach to studying memory effects is purely phenomenological: 
memory variables in the restitution relationships are usually postulated and their properties are 
obtained inductively from measurements of the differences between results from various restitution 
protocols. Asymptotic analysis such as the one used in the present work offers a deductive ab initio 
way of treatment of the memory variables. In such setting, the number of memory variables equals 
dimensionality of the phase space of the slow subsystem minus one, and the memory variables 
themselves may be identified with values of the slow variables apart from E, measured during an 
excitation front. For example, in the Caricature Noble model here there would be one memory 
variable, say taken to be n*, and in the Beeler-Reuter model, there would be dimy = 5 of them, 
and they may be identified with the values of the slow variables j, xi d, f and [Ca] as measured at 
the front, thus making the restitution much more difficult to predict. However, empirical evidence 
from simulations of modern detailed ionic models suggests that variety of slow trajectories may be 
much narrower, i.e. they may de facto be restricted to a manifold of a smaller dimensionality [55]. 
A possible theoretical explanation for such dimensionality reduction is the presence of further small 
parameters within the slow subsystem; this mechanism was considered in detail in [9] . So detailed 
studies of finer asymptotic structure of slow subsystems of practically interesting ionic models are 
a promising direction for further studies, which may help in understanding the memory effects in 
restitution. It is well known that memory effects can play considerable part in development of 
alternans and therefore in development of cardiac arrhythmias (see e.g. [43] for a recent study) 
and are for this reason of considerable interest. 

^Echebarria and Karma [21] considered a simplified model with some features of the ionic models and similar to 
the Caricature Noble model considered there, but only two variables, without any slow gates. 
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Perspective: a novel numerical method for excitation propagation The above consid- 
erations outline possible qualitative applications of the present study. A possible quantitative 
application is an advanced method for calculation of activation sequences, which can be achieved 
by the aforementioned coupling of differential equations for the local slow dynamics and for the 
front motion. For these applications, two or three spatial dimensions rather than one arc inter- 
esting, hence the equation of motion for the front is a partial differential equation of motion of a 
line (in 2D) or surface (in 3D). One immediate difference is that propagation of the front shall no 
longer depend only on the prefront voltage and j-gate, but also on the spatial configuration of the 
front. It is well known that unless the shape of the front deviates very strongly from plain, the 
effect of its shape is mostly accounted for through its mean curvature, and we have demonstrated 
that the effect of curvature can be easily incorporated into the asymptotic description of the front 
dynamics [55]. This approach can be used to describe normal activation sequences in the heart, 
when the graph of the front solution in the space-time is a manifold without internal boundaries. 
More serious problems occur if there are propagation blocks and/or wave breaks, which introduce 
boundaries of the front manifold in space-time. In such cases, a separate asymptotic description 
for the codimension-two areas, the wave break trajectories and the propagation block loci, are 
needed; obtaining such asymptotic description is another important direction for further research. 

An obvious caveat here is, as with any asymptotic approach, that the asymptotics are valid 
in the limit of e — > whereas we apply it for a finite value of that small parameter. Hence 
applicability of this approach on the quantitative level will depend on whether the error generated 
by such approximation is tolerable for the particular application; however, it should be remembered 
that higher-order terms if necessary can improve the accuracy, see an example in [9] . 
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